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During  the  second  year  of  funding  under  this  grant  substantial  progress  has  been 
made  in  the  area  of  optical  flow  computations.  Following  a  detailed  survey  of  optical 
methods  for  a  specially  designed  set  of  target  tracking  problems,  a  new  iterative  re¬ 
finement  procedure  was  developed  that  compensates  for  errors  in  the  estimated  flow 
velocity  induced  by  finite  difference  approximations  of  spatial  and  temporal  image 
derivatives.  When  coupled  with  a  total  least  squares  (TLS)  solution  scheme  (as  op¬ 
posed  to  a  more  standard  least  squares  method)  and  an  aggregate  velocity  assignment 
for  distant  targets,  this  approach  yields  very  accurate  flow  estimates  at  a  fraction  of 
the  computational  cost  of  competing  methods.  This  approach  has  the  pleasant  fea¬ 
ture  that  it  admits  a  simple  convergence  analysis  showing  that  temporal  and  spatial 
smoothing  improves  the  iterative  convergence  for  the  special  case  of  convective  flow. 
This  work  has  recently  been  supplemented  with  a  ‘‘size  band  limited”  image  process¬ 
ing  filter  based  on  partial  differential  operators;  this  approach  has  the  advantage  that 
it  is  considerably  simpler  than  a  similar  scheme  based  on  morphilogical  set  filters. 

The  publications  list  is  followed  by  a  short  narrative  section  that  describes  the 
research  that  has  been  done  so  far  under  this  grant  and  our  plans  for  the  coming  year. 
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Brief  descriptions  of  the  research  performed  under  ONR  Grant  N00014-92-J-1706  are 
given  in  this  section  with  citations  keyed  to  the  above  list.  A  more  detailed  exposition 
can  be  fotmd  in  the  original  grant  request. 
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1)  Optical  Flow  Computation:  In  [1]  the  results  are  given  of  a  long  term  project 
aimed  at  uniform  evaluation  of  optical  flow  algorithms  on  a  set  of  problems  related 
to  target  tracking  for  missiles.  The  problem  set  contains  a  variety  of  parameter 
dependent  problems  involved  in  acquisition  and  tracking  after  acquisition  especially 
in  the  presence  of  cluttered  environments.  The  set  of  basic  algorithms  tested  includes 
normal  flow  subject  to  smoothness  constraints  in  the  manner  of  Horn  and  Schunck, 
edge  dectection  via  zero  crossings  as  per  Duncan  and  Chou,  estimating  multiple 
component  motion  as  per  Bergen  et  al.  and  window  matching  across  frames  using 
differencing,  correlation  or  more  computational  intensive  distance  measures  such  as 
the  Hausdorff  distance  between  sets. 

Each  of  the  algorithms  tested  was  been  coupled  with  a  multi-resolution  scheme 
in  which  the  individual  sequence  images  are  subjected  to  a  scale  space  type  image 
processing.  The  optical  flow  is  determined  at  each  scale  with  subsequent  comparison 
of  flow  results  across  scales.  Multi-resolution  schemes  generally  fall  into  two  major 
classes  consisting  of  wavelet  decompositions  and  scale  space  decomposition  based  on 
the  evolution  of  the  image  under  a  partial  differential  operator.  A  nice  example  of  the 
latter  class  is  the  commonly  used  scale  space  resulting  from  Gaussian  smoothing  with 
the  smoothing  variance  corresponding  to  the  scale  variable.  This  in  turn  is  equivalent 
to  image  evolution  imder  the  heat  operator.  More  recently  an  highly  nonlinear  partial 
differential  operator  has  been  used  in  afl^e  invariant  image  processing.  Related 
operators  which  depend  on  the  curvature  of  the  image  isoclines  have  been  included 
in  the  scale  space  testing. 

Details  on  the  optical  flow  algorithms  and  multiscale  resolutions  are  described  in 
[1],  edong  with  numerical  resxilts  and  plans  for  future  research  in  this  area. 

Smoothness  constraints  for  optical  flow  often  involve  computations  over  a  2-D 
region,  which  may  be  the  entire  image  frame.  Hildreth  has  proposed  an  optical  flow 
computation  based  on  using  the  edges  in  an  image  sequence.  In  essence  this  is  a 
1-D  version  of  Horn  and  Schunck’s  algorithm.  Sundareswaran  and  Mallat  combine 
Hildreth’s  approach  with  the  multiscale  information  to  regularize  the  optical  flow 
computation.  Edge  based  optical  flow  methods  have  several  advantages:  (1)  they 
have  the  potential  to  overcome  the  apertme  problem  that  results  from  the  underde- 
termined  nature  of  the  constant  brightness  assumption,  and  (2)  they  greatly  reduce 
the  computational  burden  over  methods  that  rely  on  global  smoothing  constraints 
since  these  usually  lead  to  linear  systems  on  the  order  of  the  number  of  pixels  in  the 
entire  image. 

In  [2]  we  discuss  Hildreth’s  approach  and  the  Sundareswaran  and  Mallat ’s  algo¬ 
rithm  as  well  as  a  new  approach  based  on  aggregate  velocity  and  iterative  refinement. 
In  this  approach  we  assume  that  the  object  of  interest  has  edge  points  with  nearly 
equal  optical  flow  vectors.  This  is  often  the  case  for  example  with  distant  objects  in  a 
targeting  environment,  at  least  in  the  initial  stages  of  tracking.  The  aggregate  veloc¬ 
ity  asstunption  yields  an  overdetermined  system  of  equations  with  two  unknowns  (the 
X  and  y  components  of  the  group  velocity).  The  matrix  of  coefl&cients  and  the  data 
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vector  for  this  overdetermined  system  may  contain  errors  as  a  result  of  approximating 
tempor2d  and  spatial  derivatives.  Because  of  this  it  is  appropriate  to  solve  this  sys¬ 
tem  via  the  total  least  squares  approach  rather  than  the  more  common  least  squares 
method.  This  is  described  in  more  detail  in  [2].  Because  there  are  only  two  unknowns, 
the  computation  of  the  aggregate  velocity  vector  is  very  fast.  The  computed  velocity 
may  still  be  in  error  because  of  the  derivative  approximation  errors.  Because  of  this 
we  introduce  an  iterative  correction  procedure  that  subtracts  the  computed  velocity 
from  the  image  sequence  and  then  recomputes  a  new  aggregate  velocity  that  in  turn 
may  be  subtracted  from  the  new  image  sequence.  This  is  repeated  until  the  velocity 
stabilizes. 

An  important  question  is  whether  indeed  this  method  does  converge.  This  problem 
admits  analysis  for  the  case  of  one  dimensional  convected  flow  and  conditions  ensur¬ 
ing  convergence  are  easily  formulated  (see  [2]).  Moreover  we  show  that  smoothing 
(either  temporal  or  spatial)  enhances  convergence  without  changing  the  flow  velocity. 
Numerical  examples  are  presented  illustrating  the  main  points  of  our  procedure  for 
both  one  and  two  dimensional  problems. 

2)  Sensitivity  and  Condition  Estimation:  One  of  the  main  goals  of  this  grant 
has  been  to  develop  a  general  procedure  for  measuring  the  sensitivity  of  wavelet  and 
signal  processing  algorithms.  Such  a  procedure  should  be  able  to  provide  an  efficient 
assessment  of  the  reliability  of  computed  values  from  these  algorithms.  Fortunately, 
the  work  in  [3]  forms  a  basis  for  this  t]rpe  of  estimation  procedure  for  multivariable 
scalar  valued  Actions  f  :  R*  R.  For  these  functions,  the  norm  of  the  gradient 
provides  a  measure  of  the  sensitivity  of  the  function  at  the  point  of  evaluation.  The 
main  problem  then  is  to  get  an  estimate  of  the  norm  of  the  gradient  without  having 
to  evaluate  the  gradient  itself,  since  this  would  require  n  function  evaluations  if  we 
use  finite  differences  to  approximate  the  n  partial  derivatives  in  the  gradient.  The 
saving  grace  of  this  problem  is  that,  by  using  the  theory  of  random  iimer  products  on 
the  unit  sphere  in  n  dimensions,  we  can  estimate  the  norm  of  the  gradient  with  only  a 
few  function  evaluations  (typically  two  or  three)  and  that  exact  probabilistic  bounds 
can  be  given  on  the  accuracy  of  the  estimate.  For  example,  condition  estimates  are 
considered  acceptable  if  they  are  within  a  factor  of  10  (because  we  are  interested  in 
estimating  the  number  of  digits  of  accuracy).  If  three  function  evaluations  are  used 
in  estimating  the  condition  number  of  a  scalar  valued  function,  then  the  probability 
that  the  estimate  is  within  a  factor  of  10  is  0.9989  and  k  function  evaluations  lead 
to  order  condition  estimates,  i.e.  the  probability  of  being  off  by  a  factor  u)  is  less 
than  a  constant  times  1/a;*.  A  very  readable  accoimt  of  this  work  has  appeared  in 

[4]. 

This  ‘‘small-sample”  statistical  condition  estimation  method  has  been  extended  to 
matrix  valued  functions  by  the  work  in  [5].  Here  there  are  two  possible  approaches: 
1)  treat  the  individual  entries  of  the  function  as  if  they  were  scalar  valued  functions 
and  apply  the  scalar  condition  estimates,  or  2)  estimate  the  overall  condition  of  the 
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functional  map  via  the  norm  of  the  Frechet  derivative.  The  first  approach  has  the 
nice  feature  that  the  scalar  theory  can  be  applied  immediately  and  that  estimates 
of  the  conditioning  of  all  the  entries  can  be  obtained  at  the  same  time  with  only  a 
few  function  evaluations.  However,  if  the  number  of  outputs  is  large  then  it  is  often 
more  convenient  to  provide  just  one  number  as  an  overall  condition  measure,  as  in 
the  second  approach.  This  overall  assessment  has  been  considered  in  [5]  in  which  a 
conjecture  is  presented  relating  the  theory  of  the  accuracy  of  the  individual  estimates 
with  that  of  the  overall  condition  estimate.  A  conservative  form  of  this  conjecture  is 
proved  in  [5],  but  the  accuracy  estimates  of  this  conservative  appraoch  are  too  pes¬ 
simistic  in  the  case  of  a  very  large  number  of  inputs.  Although  extensive  numerical 
tests  show  that  it  applies  in  many  example  problems,  more  work  is  needed  in  this 
area  to  establish  the  full  conjecture. 


3)  Wavelets  and  Differential  Equations:  Wavelets  and  differential  equations 
have  natural  connections  from  two  points  of  view:  1)  many  physical  processes  have 
self  similar  properties  under  translation  and  rescaling,  and  2)  recent  work  (L.  Al¬ 
varez,  F.  Guichard,  P.  Lions  and  J.  Morel,  ‘‘Axiomatic  and  Fundamental  Equations 
of  Image  Processing,”  preprint)  shows  that  the  image  processing  schemes  which  are 
invariant  with  respect  to  certain  image  operations,  such  as  transposition  or  rotation, 
can  be  viewed  as  time  slice  projections  of  differential  operators.  As  an  example,  the 
heat  operator  gives  rise  to  smoothing  or  blurring  image  processors  that  are  equivalent 
to  Gaussian  filtering;  Some  care  must  be  taken  in  the  numerical  solution  of  affine 
invariant  nonlinear  differential  equations;  the  work  in  [6]  shows  that  standard  finite 
difference  implementations  can  produce  erroneous  results  near  extrema.  The  family 
of  exact  solutions  in  [6]  points  the  way  to  avoiding  these  problems  ^d  at  the  same 
time  leads  to  a  nonlinear  version  of  the  Courant  stability  relation  between  the  spatial 
grid  size  and  the  scaling  step  size. 

The  affine  partial  differential  operator  appears  to  be  particularly  useful  in  filtering 
that  is  band  limited  with  regard  to  object  size  using  tophat  operators  based  on  PDEs. 
In  this  procedure,  evolution  of  the  image  tmder  the  PDE  operator  is  identified  with 
morphological  erosion  and  backwards  time  evolution  is  identified  with  dilation.  The 
composition  of  these  operations  is  subtracted  from  the  original  image  to  give  the 
tophat  operator.  See  [7]  for  a  discussion  of  tophat  operators  and  morphological  filters. 

More  specifically,  let  u  —  u{x,  y,  0)  be  an  image  and  define  the  scale  space  resolu¬ 
tion  of  u  as  the  solution  u  s  u(z,  y,  s)  to  the  initial  value  PDE 


du 

dt 


=  Lu, 


(1) 


with  u{x,y)  at  s  =  0  equal  to  u(3;,  y,0).  Here  L  is  a  spatial  partial  differential 
operator.  For  the  affine  invariant  operator 

•LaAne  ~  ”  2u*yU*tZy)^^^.  (2) 
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Let  F,{w)  denote  an  image  w  after  it  has  evolved  for  a  time  period  s  under  the 
PDE.  The  tophat  operator  is  defined  by 

r.(u)  =  u  -  (3) 

That  is  first  we  let  u  evolve  for  time  s,  then  we  reverse  the  evolution  and  subtract 
this  result  from  the  original  image. 

To  provide  some  motivation,  we  note  that  under  the  affine  invariant  tophat,  a 
peak  of  radius  r  with  circular  symmetry  will  disappear  at  s  =  since  this  is  how 
long  it  takes  the  circular  level  curve  of  radius  r  to  move  inward  to  the  center  of  the 
peak.  Thus  any  isolated  peak  in  the  original  image  that  is  of  radius  less  than  or  equal 
to  r  will  disappear  by  the  time  s  =  When  the  affine  PDE  is  then  reversed  and 
run  backward  small  peaks  (of  radius  <  r)  are  not  recreated  because  the  information 
necessary  to  do  the  restoration  has  disappeared.  Thus  these  peaks  are  missing  from 
F-,{F,{u)),  while  peaks  of  larger  radii  are  restored.  Note:  only  the  outer  rings  of 
these  larger  pe2dcs  are  restored;  the  inner  parts  of  these  peaks  are  destroyed  just  as 
the  smaller  peaks  are.  In  particular  the  tophat  operator  uses  this  selective  destruction 
to  eliminate  peaks  beyond  a  given  radius;  by  differencing  tophat  operators  of  different 
radii  we  obtain  size  limited  band  filters. 

The  connection  between  image  processing  filters  and  evolution  under  a  partial  dif¬ 
ferential  operator  opens  up  the  possibility  of  studying  the  sensitivity  of  wavelet  and 
signal  processing  algorithms  via  the  sensitivity  of  solutions  to  differential  equations 
and  vice  versa.  The  general  question  of  the  sensitivity  of  solutions  to  differential 
equations  has  also  taken  on  increased  importance  to  researchers  at  the  NAWC  be¬ 
cause  of  the  problem  of  obtaining  accurate  computer  predictions  of  variables  related 
to  radome  electromagnetic  problems,  such  as  the  boresight  error  problem.  We  are 
especially  interested  in  extending  the  sensitivity  analysis  of  [3]- [5]  to  the  methods 
developed  in  [8]- [10]  for  solving  a  variety  of  problems  in  differential  equations.  Like¬ 
wise  the  work  in  [11]  opens  up  new  avenues  for  sensitivity  estimation  in  differential 
problems  related  to  control  theory. 

4)  Research  Directions  for  1994-1995:  Optical  Flow  Methods  Optical 
flow  algorithms  for  detecting  motion  in  cluttered  environments  have  assumed  new 
importance  in  missile  tracking  problems  at  the  Naval  Air  Warfare  Center  in  China 
Lake,  especially  in  the  area  of  thermal  imaging  against  a  hot  background  as  opposed 
to  ‘blue-sky’  tracking.  The  essential  problem  to  be  treated  is  the  reduction  in  the 
amount  of  computation  needed  to  obtain  smooth  flow  fields  via  global  constraints 
on  the  allowed  variation  in  the  field.  Work  on  this  problem  has  shown  that  the 
iterative  smoothing  procedures  needed  to  solve  the  globad  coupling  problems  can  be 
greately  speeded  up  by  working  in  a  cascade  of  resolution  levels  such  as  those  provided 
naturally  by  wavelet  decompositions. 

A  second  approach  to  the  problem  of  computational  reduction  involves  translating 
the  two-dimensional  imaging  problem  into  a  one-dimensional  problem  by  focusing 
only  on  image  aspects  such  as  edges  and  comers.  Thus  the  identification  of  such 
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structures  for  example  by  tophat  operators  takes  on  new  urgency.  An  important 
problem  in  this  area  is  the  invariance  of  the  computed  optical  flow  field  with  respect 
to  reflections  in  the  x  and  y  axes;  invariance  of  this  type  seems  to  depend  on  the  choice 
of  the  wavelet  filters  that  are  used  in  the  processing  especially  with  regard  to  their 
symmetry  with  respect  to  reflections  in  x  and  y.  Research  in  this  area  will  focus  on 
eliminating  invariance  problems  from  existing  edge  detection  optical  flow  algorithms 
and  enhancing  their  resolution  by  modifying  the  wavelet  filters  to  correspond  to  the 
image,  i.e.  using  the  techniques  of  selecting  the  best  wavelet  to  improve  the  optical 
flow  calculations. 

Research  in  this  area  will  be  conducted  in  conjimction  with  an  ongoing  effort  at 
NAWC  under  the  direction  of  Dr.  Gary  Hewer. 


Abstracts 

The  following  pages  give  abstracts  from  papers  accepted  or  submitted  under  ONR 
Grant  Ntunber  N00014>92*J-1706. 
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A  Survey  of  Optical  Flow  Methods  for  TVacking  Problems 

Charles  Kenney,  Gary  Hewer  and  Wei  Kuo 
Naval  Air  Warfare  Center 
China  Lake  CA  93555-6001 

ABSTRAdT 

Itoulta  ara  predated  of  a  aumetieal  ,a„ay  of  optical  Jo-  algorithm,  fo,  ttacktag  probUm.  a..ociat«i  with 
inftated  .magingin  high  .peed  mWM.  The  algorithm,  tmtad  mclud.  th<».  aa-^iatad  with  normal  dow  .ahject 
to  global  .raoothne»  conatrainu,  adg.  detection  via  .ero  ctoMlng.  of  the  image  after  convolution  with  .patio- 
t«nporai  Jlter.,  and  windowed  matching  technique..  The  tracking  problem,  conridered  in  thi.  .urvey  fall  into  two 
cl»,.,,  acquMtion  and  tracking  after  acquieition.  There  cl^we.  can  be  further  divided  into  near  and  far  range 
problem.,  charac.eri.ed  by  extended  and  point  target  image..  Other  parmneter.  of  intereet  model  allowable  target 
and  .ensor  motion  and  the  amount  of  background  clutter.  Ewi  of  the  above  method,  fo.  determining  optical  do- 
can  be  need  in  conjunction  with  a  variety  of  image  p..proc.»ing  technique,  .uch  a.  kernd  wnoothing  (eepecially 
by  Gau^ian  kerneb),  and  evolution  under  afflne  invariant  partial  differ^dial  operator..  Th.«  preproewwing 
method,  cm.  J,o  be  u«d  in  combination  with  other  approach..  «.d.  a.  temporal  layering  in  which  euccive 
.mage,  are  comb.nml  to  produce  tamje.  with  .treak.  who.,  edge,  are  predominantly  paralW  to  the  optical  dow. 

2.  INTRODTinTTON 

In  thi,  paper  w.  report  on  a  long  term  projw:t  aimed  a.  uniform  evaluation  of  optical  dow  algorithm,  on 
a  «t  of  problem,  relatml  m  ..,g„  ...dring  fo,  mimile,.  The  pneblmn  «.  oontain.  a  variety  of  parameter 
ependm..  problem,  .nvolved  in  «qui.i.ion  and  tracking  after  acquirition  eepecially  in  the  pr.,«.c.  of  clu.teed 
euvnonmenu.  The  .e.  of  haaie  algorithm.  teW  „  fa,  include,  normal  do-  object  ri,  mioothnem  conatraint,  in 
t  e  manner  of  Horn  and  Shunk  (15),  edge  dectectiou  via  nro  erouing.  a.  per  Duncan  and  Chou  [9],  eetimating 
multiple  component  motion  a.  per  Bergmi  e.  al.  [g],  and  window  matching  acorn,  ftame.  uring  differencing, 
»rr.lat.on  or  more  computation^  intenriv.  diatanc.  memure.  .uch  a.  the  Hauadorff  diatmice  betwemi 
6].  The  «t  of  opticj  dow  algoriUnn,  will  he  mtpmided  in  the  future  to  include  other  approtohe.  ,uch  the 
^neraliaation,  by  Schnor,  [30]  of  toe  „oo.hn,»  funCionab  of  Horn  and  Schunck  (Id)  and  Nagri  and  Enkelman. 
(23),  Kalman  dltormg  .imwUecollfaione,.imation  of  [22]  and  other  Kalman  ptocedure,  (31),  «lg,  detection  baaed 
on  phe«  mformation  (11)  a.  weU  »  aldne  invariant  movie  mialyri.  method,  (12). 

Each  of  the  algorithm,  toatod  ha.  been  coupW  with  a  multi-rewiluaon  tohmne  in  which  the  individual 
toquence  miage.  are  .ubjected  to  a  male  .pace  type  image  procecing.  The  opticri  dow  b  determined  at  each 
.cal.  with  „bto,uen.  comparbon  of  dow  reeulb  acme,  aerie,.  Multi-rewilution  toheme.  gmierrily  fril  into  two 
major  dame,  conamring  of  wavelet  decompoaition,  mid  aerie  .pm,  deeomporition  haaml  on  the  evolution  of  the 
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ABSTRACT 


Optical  flow  is  an  estimate  of  the  velocity  fleld  based  on  the  change  of  intensity  patterns 
in  successive  images,  and  is  an  important  quantity  in  computational  vision  for  dense  images. 
Because  of  the  aperture  problem  optical  flow  computations  can  be  ill-posed.  This  problem 
is  compounded  by  derivative  estimation  errors.  This  paper  presents  an  aggregate  velocity 
scheme  that  uses  iterative  velocity  refinement  along  object  edge  contours  obtained  via  the 
Mallat-Zhong-Hwang  wavelet  and  chaining  algorithms.  By  working  with  edge  information 
and  aggregate  velocities  we  avoid  the  aperture  problem;  iterative  refinement  compensates 
for  errors  in  the  derivative  estimation.  Our  approach  assigns  a  common  velocity  to  the 
edge  points  of  an  image.  When  combined  with  a  constant  brightness  assumption  this  yields 
an  overdetermined  set  of  linear  equations.  Since  the  data  vector  and  matrix  coefficients 
of  this  linear  system  consist  of  temporal  and  spatial  derivative  estimates,  respectively,  and 
both  are  subject  to  errors,  the  overdetermined  system  is  solved  using  a  total  least  squares 
approach.  The  resulting  velocity  estimate  is  then  subtracted  from  the  image  sequence  and 
the  velocity  estimation  procedure  is  repeated  for  the  new  image  sequence.  This  approach 
is  very  fast  and  accurate  for  images  that  have  nearly  the  same  edge  velocity  vectors  as  is 
usually  the  case  for  distant  objects.  A  convergence  analysis  is  given  for  the  special  case 
of  one  dimensional  convected  flow  and  it  is  shown  that  spatial  and/or  temporal  smoothing 
enhances  the  convergence. 


I.  INTRODUCTION 


Optical  flow  is  an  estimate  of  the  velocity  field  based  on  the  change  of  intensity  patterns 
in  successive  images,  and  is  an  important  quantity  in  computational  vision  in  dense  images. 
A  well-known  method  for  optical  flow  computation  is  the  optimmation  approach  introduced 
by  Horn  and  Schunck  [4].  In  this  approach,  the  sum  of  two  cost  functionals  of  the  unknown 
flow  vectors  is  minimized.  One  of  the  cost  functional  is  the  ‘‘brightness  constraint”.  The 

’Approved  for  public  release;  distribution  is  unlimited. 
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Small-Sample  Monte  Carlo  Condition  Estimates 
for  General  Matrix  Functions  * 
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Abstract 

A  new  condition  estimation  procedure  for  general  matrix  functions  is  presented  that  ac¬ 
curately  gauges  sensitivity  by  measuring  the  effect  of  random  perturbations  at  the  point  of 
evaluation.  The  number  of  extra  function  evaluations  used  to  evaluate  the  condition  estimate 
determines  the  order  of  the  estimate.  That  is,  the  probability  that  the  estimate  is  off  by  a  given 
futor  is  inversely  proportional  to  the  factor  raised  to  the  order  of  the  method.  Comparison  is 
given  with  the  power  method  of  condition  estimation  which  because  of  its  sequential  nature 
potentially  require  twice  the  computational  effort,  as  is  the  case  for  the  problem  of  estimating 
the  senntiviy  of  the  matrix  exponential.  A  group  of  examplm  illustrates  the  flexibility  of  the 
new  estimation  procedure  in  handling  a  variety  of  problems  and  tjrpes  of  sensitivity  estimates, 
such  as  mixed  and  componentwise  condition  estimates. 

Keifwordr.  Monte  Carlo,  conditioning,  matrix  functions 
AMS(MOS)  nbject  cla$$ifieaiion:  6SF35,  05F30, 1SA12 
Ahhnviaiei  titk:  Monte  Carlo  Condition  Estimates 


1  Introduction 

If  a  scalar  real-valued  function  /  has  a  large  derivative  at  x  €  11  then  the  rdationship 

/(x  +  «)*/(x)  +  /'(x)fl  +  0(f») 

shows  that  small  perturbations  in  x  can  lead  to  relatively  large  changes  in  /.  The  idea  that  the 
magnitude  of  the  derivative  provides  a  measure  of  the  sensitivity  or  conditioning  of  /  carries  over 
to  higher-dimensional  mappings  [5],  [21],  [31],  [29],  [32],  [28],  [35],  [40].  For  example,  if  /  maps 
E,"  into  E,  the  gradient  of  /  at  x  is  the  row  vector  V/(x)  *  (6/(x)/dxi, . .  .,d/(x)/5xO  and  the 
Taylor  expansion  of  /  has  the  form 

_ f{t  +  S*)  =  /(»)  +  SVfix)z  +  0(0  .  (1) 

'This  Nseuch  wts  supported  in  pert  by  the  N&tioul  Sdeace  Fonadutiou  (tad  AFOSR)  under  Grunt  No.  ECS87- 
18897,  the  Air  Force  Office  of  Sdeatffic  Research  under  Contract  No.  AFOSR-91-0340,  and  the  Office  of  Naval 
Research  under  Grant  No.  N00014-93-J-1706. 
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Abtireet 

Ibk  pap«  •  nent  itetiftieally  baitd  mtae-aenn  «ftiiiiator  to  metriMi.  Tb*  a«w 

timetof  leqoini  oaljr  e  ftw  uutrix- vtetor  multiplicetioM  aad  caa  b«  appUad  wbae  tba  matrix 
ia  net  known  c^Udtiy.  It  ia  naaAil  tot  afldantbr  artfanating  the  amrithritjr  ot  vaetoe-valnad 
flmctioai  aad  can  be  applied  to  many  pioblema  wham  the  power  method  mae  into  difflcattka. 
Low«  bonade  br  the  piebabOity  that  aa  ertimate  ia  within  a  gim  (hetor  ci  the  eocieet  norm 
ate  derived.  Three  hmda  are  atrairittferward  to  eempnte  aad  ahew  that  a  vaqr  iaaeenrate 
aatimate  ia  aactiamaly  naUhafy  ia  moat  eaeaa.  A  eooaarvative  low«  bovad  hae  bean  derived  aad 
a  tiriiter  bonad  ia  given  ia  tte  fam  of  a  conjeetnia.  Ihia  coe|ectaie  ia  tina  ia  aonm  iB^ortaat 
ape^  eaeaa  aad  the  gaaarai  caee  »  at^pertod  by  cooaidarabk  oat^tUai  evidanee. 

ifaywerda:  Cenditjonfaig,  amtrix  Aiaetioaa 
AMS(MOS)  toijtet  oktoi/UtHoK  dSFSS,  esrSO.  1U19 
AUnoiotoi  HUm  Statiatieal  Coaditian  Eatimataa 


1  Introduction 

A  novel  method  for  effldeatly  eatimatiBg  the  aenaitivity  of  a  acaler  Ibaetion  at  a  point  waa  recentfy 
intxodncod  in  (17).  Thla  nMthod  la  baaed  on  Implicitly  piojecring  an  appradmate  gradient  of  the 
function  onto  a  uaBbnnly  randomly  doaea  low-order  anbapaeo,  aad  than  compntiag  the  norm  of 
the  remit.  Properly  acalod,  thla  givea  an  eatfanate  of  the  aoim  of  tho  gradient,  aad  thereby  the 
condition  nnaher  of  the  fhnetian  at  the  eatimarion  point.  The  method  reqataea  only  tho  evalnatkm 
of  the  Ibaetion  at  thia  point  aad  a  bw  nearby  pointa. 

Thla  p^er  aactenda  the  reanha  (tf  [17]  to  the  eatimation  of  the  Ibobeahu  norm  of  the  Jacobian 
of  a  geaanl  veelot-vntnod  fbnction.  Thna,  let  / :  K*  E*  be  a  fhnction  dUbrentUMa  at  a  point 
X,  aad  define  the  JneoUan  at  s  aa 

*TUeNMaMhwwMppoctadbytheN«tiaaalSdeMaIbeadatkeeadar<3caatlf»  ECS4i]0a4S,tkeAirl<Dtev 
Oflka  of  Sdamile  Raaaaidi  aadw  Coattaet  No.  ATOSlUai-Otao,  aad  the  Oftee  af  Naval  Baaaaich  eadar  Ceatnct 
No.  N00014-f3>MTaa. 
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Abstract 

A  new  approach  is  presented  for  estimating  the  conditioning  of  general  matrix 
functions  by  measuring  the  effect  of  random  perturbations  at  the  point  of 
evaluation.  This  method  is  efficient  in  the  sense  that  the  number  of  extra 
function  evaluations  used  to  evaluate  the  condition  estimate  determines  the 
order  of  the  estimate.  That  is,  the  probability  that  the  estimate  is  off  by 
a  given  factor  is  inversely  proportional  to  the  factor  raised  to  the  order  of 
the  method.  The  “transpose-free”  nature  of  this  new  method  allows  it  to  be 
applied  to  a  much  broader  range  of  problems  than  the  commonly  used  power 
method  of  condition  estimation.  A  group  of  examples  illustrates  the  flexibility 
of  the  new  estimation  procedure  in  handling  a.variety  of  problems  and  types  of 
sensitivity  estimates,  such  as  mixed  and  componentwise  condition  estimates. 
Short  Matlab  routines  are  included  to  demons^ate  the  ease  with  which  the 
new  condition  method  can  be  implemented  in  a  general  setting. 


1 .  Introduction 

If  a  scalar  real-valued  function  /  has  a  large  derivative  at  r  €  IR,  then  small 
perturbations  in  x  lead  to  relatively  Urge  changes  in  /,  The  idea  that  the 
magnitude  of  the  derivative  provides  a  measure  of  the  sensitivity  or  condition* 
ing  of  /  carries  over  to  higher-dimensional  maps  [3],  (14),  [22],  [20],  [23],  [19], 
[28],  [31].  For  example,  if  /  maps  IR"  into  R,  the  gradient  of  /  at  *  6  IR"  is 
the  row  vector  7/(x)  =  (d/(x)/dxi . a/(x)/a*n).  and 

/(i  Sz)  =  /(x)  +  6Vf(x)z  +  O(S^)  .  (1) 

where  ]|r||  =  i  with  ||  •  ||  denoting  the  vector  2-norm  defined  by  ||z||’  =  E|s,p, 
and  j  is  a  small  positive  real  number.  As  a  matter  of  notation,  we  use  to 
denote  the  gradient  ^ 

»^  =  y/(z).  (2) 

It  can  be  seen  from  (1)  that  the  norm  of  v  is  an  appropriate  measure  of 
tlie  local  sensitivity  of  /.  For  simplicity  of  exposition,  we  assume  through¬ 
out  this  paper  that  all  functions  considered  are  at  least  twice  continuously 
differentiable. 

Somewhat  surprisingly,  rigorous  probability  arguments  show  that  only  a 
few  function  evaluations  are  needed  to  obtain  accurate  and  reliable  estimates 
of  the  norm  of  v.  In  this  paper  the  small-sample  method  of  [24]  is  applied  to 
the  problem  of  condition  estimation  for  more  general  matrix  functions  such 
as  the  matrix  exponential  map  X  e^*^  and  the  map  (A,  F,C)  X.  where 
X  satisfies  an  algebraic  Rtccati  equation 

0  =  C7+  A^X  -h  XA  -  XrX.  (3) 


This  is  illustrated  in  Section  2  with  many  other  examples.  That  section  also 
compares  the  small-sampU statistical  method  with  the  more  traditional  power 
method  of  condition  estimation.  Short  Matlas  routines  are  given  in  Section 
2  for  estimating  the  sensitivity  oC  general  matrix  fiiactione.  The  theoretical 
basis  for  the  small-sample  methods  is  briefly  described  in  the  remainder  of 
this  section  and  is  given  in  more  detail  in  [24]. 

The  central  idea  of  this  method  is  that  the  norm  of  the  gradient  can  be 
estimated  if  we  can  afford  another  function  evaluation  beyond  /(x).  Suppose, 
for  example,  that  we  evaluate  /  at  x+6z.  If  z  has  unit  norm,  then  the  fiewton 
difference  defined  by 


dx  = 


/(x-hSz)^f(x) 

S 


(4) 


'This  res«u^  wm  supporud  in  put  by  the  Nniionri  Sdcnc*  Foundation  (rant  ECS- 
9120643,  th«  Air  Force  Office  of  Scientific  Reeearch  (rant  AFOSR-9 1*0340,  the  Office 
of  Navel  jRcsearch  grant  N00014>93-J-1706. 


satisfies 

dz  =:  v^x  -1-  0(6). 

(5) 

If  X  is  selected  uniformly  and  randomly  from  the  unit  sphere  5„_i 
mensions,  then  (see  Theorem  1  in  [24])  the  expected  value  of  |u^;| 
by 

in  n  di- 
is  given 

=  ||v||X„ 

(6) 

where  S,  =  1,  and  for  n  >  1 

(T) 

p  _  2  2.4.6  -  (n-2)  , 

Z  1  O  e  /  lor  n  even. 

V  1  -3  •5-“(n  -  1) 

(3) 

For  large  n  it  is  convenient  to  use  the  approximation  £«  ss  v/2/(T(n  -  1/2)). 
which  is  accurate  to  three  significant  digiu  for  n  >  10. 

The  condition  estimator  |dx|/£n  has  expected  value  equal  to  the  true  con¬ 
dition  number  ||v||  plus  a  term  of  order  S.  Typically,  we  can  take  6  sufficiently 
small  (say,  less  than  lO’*).  to  that  for  the  purposes  of  estimating  ||i;i|  the  ex¬ 
pected  value  of  \dx\/ En  is  more  than  adequate.  To  avoid  having  to  add  terms 
of  order  ^  to  all  of  our  results,  we  will  now  make  the  simplifying  assumption 
that  for  any  given  vector  x  €  we  are  able  to  evaluate  w^x  with  the 
understanding  that  in  real  life  we  are  merely  able  to  evaluate  dx  as  in  (4). 
This  is  reasonable  when  we  remember  that  for  most  sensitivity  estimates,  we 
only  need  to  know  the  true  condition  number  to  within  a  factor  of  10  or  so, 
and  we  can  often  tolerate  errors  in  the  estimate  up  to  a  factor  of  100. 

What  is  the  probability  that  the  estimator 


c  s  |v^x|/£, 

lies  within  a  given  factor  w  of  ||v||?  This  question  was  considered  in  detail  in 
[24]  and  it  was  found  that  because  has  a  Beta  distribution. 


Pr(||vl|/w  <  <  <  u;||v||)  >  1  -  -L  +  o(-^r)  .  (9) 

TU/  Wr 

Thus  (  is  a  linear  or  first-order  condition  estimate  in  the  sense  that  the  chance 
of  a  catastrophically  low  or  high  estimate  is  inversely  proportional  to  the  size 
of  the  error.  For  example,  Pr(||v||/100  <  C  <  100||v||)  >  0.9938.  so  that  the 
chance  of  being  off  by  more  tbiut  a  factor  of  100  is  less  than  1  in  100. 

While  this  is  good,  there  are  some  situations  in  which  we  need  more  reli¬ 
ability.  One  way  to  achieve  this  is  to  use  more  function  evaluations  to  get 
different  value*  corresponding  to  independently  randomly 

generated  veeton  i(‘),  r<»), . . . ,  x<*)  in  5n.,  and  then  to  take  the  average 


<(m)5 


m 


(10) 


For  V  >  1,  we  have  from  [24]  that 


Pr(IHI/w  <  <(m)  <  HMD  >  1  -  A  f— )"  +  O(Arr)  .  (11) 

m!  \TU>y 

with  asymptotic  equality  as  w  —  +oo  or  m  —  +oo.  Thus  <(m)  is  an  mth- 
order  condition  estimator.  That  is,  the  probabUity  that  the  estimator  is  off 
by  more  than  a  factor  of  w  is  leM  than  a  constant  times  w""".  For  example, 
with  m  s  3,  which  corresponds  to  three  extra  function  evaduations,  we  have 
PKiWII/10  <  C(3)  <  10||v||)  >  0.99884,  so  that  the  chance  of  being  off  by 
a  factor  of  10  or  more  is  about  1  in  862.  Similarly,  the  chance  of  being  off 
by  a  factor  of  100  or  more  is  approximately  one  in  a  million.  Hence  very 
reliable  condition  estimates  can  be  obtained  with  only  a  few  extra  function 
evaluations. 
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Abstract 

Exact  solutions  to  a  FDE  associated  with  afiSne  invariant  image  processing  are 
derived  for  initial  conditions  corresponding  to  extrema  and  saddle  points.  For  these 
initial  conditions  the  original  nonlinear  PDE  reduces  to  the  standard  convection  equa¬ 
tion.  Application  of  the  affine  invariance  property  then  generates  a  much  larg»  family 
of  exact  solutions,  including  those  with  elliptically  symmetric  initial  conditions.  Ex¬ 
amination  of  this  family  clarifies  surface  evolution  near  extrema  and  saddle  points  in 
terms  of  flows  along  characteristics:  surface  information  flows  toward  extrema  and 
away  from  saddle  points.  This  viewpoint  also  leads  to  a  nonlinear  Courant  number 
relating  the  permissable  evolution  step  size  to  a  ^ven  spatial  step  size.  The  exact  so¬ 
lution  for  circularly  symmetric  initial  conditions  shows  that  standard  central  difference 
approximations  of  the  affine  invariant  PDE  can  produce  incorrect  results  in  the  form 
of  spikes  at  extrema.  A  procedure  for  avoiding  this  problem  is  discussed  and  short 
Matlab  routines  for  affine  image  processing  are  pven. 

‘This  research  was  supported  by  the  Office  of  Naval  Research  under  ONR  Grant  Number  N00014-92-J- 
1706. 
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Abstract 

A  simple  procedure  is  given  for  constructing  a  complete  basis  of  solutions 
to  a  given  linear  constant  coefficient  partial  differential  equation.  A  variety  of 
examples  are  presented  including  eigenvalues  problems  such  as  the  Helmholtz 
equation. 


*Thia  research  was  supported  by  the  Office  of  Naval  Research  under  ONR  Grant  Number  N00014>92-J- 
1706. 
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Abstract 

Numerical  methods  are  presented  for  handling  two  of  the  most  important  computa¬ 
tional  problenu  associated  with  the  study  of  parameter  dependence  solutions  to  the 
periodically  time  dependent  Schrodinger  equation  as  approximated  by  a  finite  system  of 
ordinary  differential  equations  (ODEs):  1)  lengthy  integration  times  due  to  the  growth 
in  parameters  as  the  number  of  states  increase  and  2)  repeated  integrations  over  a  grid 
in  parameter  space.  These  methods  are  designed  to  work  together  and  are  independent 
of  the  numerical  algorithm  used  to  integrate  the  ODE  system.  The  first  method  uses 
a  simple  diagonal  transformation  of  the  solution  to  contnd  the  magnitude  of  terms  ap¬ 
pearing  in  the  differential  system,  and  the  second  method  r^es  on  matrix  interpolation 
together  with  a  Thgdor's  series  expansion  to  reduce  the  number  of  integrations  in  pa¬ 
rameter  space.  The  comUnation  of  these  two  procedures  results  in  a  dramatic  decrease 
in  computatiott  time  for  an  example  problem  of  an  riectron  confined  to  a  quantum  well. 

Keytpordr.  Schrddlnger’s  Equation,  Parameter  Study,  Matrix  Interpolation 
AMS(MOS)  iuhjtet  eUunfieationx  65L05, 65F30. 
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*Tlus  icteuA  wu  rappoitcd  by  s  Fsodoi-Lynca  RswsiA  grut  from  tlw  Akxaadct  voa  HumboUt  -  SUftvng. 

'This  nteudi  wts  snppoitsd  ia  put  bjr  the  OiSc*  of  NstsI  Rsacsidi  udu  ONK  Qnuit  Nuabu  N00014-93-J-1708. 


1 


Cncum  SYsma  Signal  Pkgcbsb 
VbL.  13.  Na  6. 199<.  Pr.  695-732 


Methods  for  the  Numerical 
Integration  of 
Hamiltonian  Systems* 

y.  J.  Hench,^  C  5.  Kenney^  and  A.  J.  Laub^ 


Abstract  To  compute  an  infinite  horizon  optimal  controller  for  a  linear  periodic  system 
via  an  invariant  subspace  method,  the  computation  of  the  period  map  associated  with  the 
Hamilton- Jacobi-Bellman  equations  is  required.  In  this  paper  we  discuss  methods  for  the 
numerical  integration  of  such  Hamiltonian  systems.  Two  numerical  integration  techniques 
are  introduced.  A  method  is  developed  whereby  sympicctic  invariants  associated  with  the 
Hamilton-Jacobi-Bellman  equations  are  preserved.  Also,  a  shifting  scheme  is  introduced 
that  in  effect  swaps  the  roles  of  the  stable  and  unstable  invariant  subspaccs  by  using  the 
semigroup  property  of  sute  transition  matrices.  A  shift  is  introduced  into  the  resultant  initial 
value  problem  ensuring  that  the  eigenvalues  of  a  differential  equation  reside  in  the  region 
of  absolute  stability  for  an  appropriate  numerical  integration  routine.  These  techniques  are 
then  compared  to  standard  numerical  integration  routines  to  ascertain  their  efficiency  and 
accuracy. 


1.  Introduction 

In  this  paper,  we  examine  various  methods  of  computing  the  optimal  feedback 
controller  gains  for  linear,  periodic  time- varying  systems 

i(0  =  Ait)x(t)  +  B(r)M(0  (1) 

with  quadratic  cost  functional: 

mini/*  [x^  Q(t)x uT  R(t)u]dt  (2) 

“2  70 
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Abstract 

Certain  wpects  of  stable  Lyapunov  operators  can  be  easily  studied  by  exploiting 
the  linearity  of  the  trace  operator  and  its  invariance  under  reversal  of  order  in  matrix 
products.  For  example,  sharp  upper  and  lower  bounds  on  the  trace  of  solutions  to  the 
stable  Lyapunov  equation  can  be  obtained  by  applying  the  trace  operator  to  a  well 
known  integral  representation  of  these  solutions. 

Other  applications  include  using  the  connection  between  dual  norms  and  the  trace 
operator  to  obtun  new  results  on  the  norms  of  Lyapunov  operators  assodated  with  the 
conditioning  of  solutions  to  the  Riccati  equation.  In  this  r^ard,  trace  norm  results  can 
be  obtained  from  wdl  known  spectral  norm  results  since  the  trace  and  spectral  norms 
are  dual  to  each  other.  A  somewhat  deeper  analypis  involving  the  power  method  gives 
monotonicaUy  decreasinf  upper  bounds  on  the  IVoboJus  norms  of  these  Ljrapunov 
operators;  these  upper  bounds  complement  the  usual  monotonicaUy  increasing  lower 
bounds  associated  with  the  power  method  and  provide  a  nice  means  of  assessing  the 
accnraqr.«f  the  resuMhf  Ikobenius  norm  estimates. 

JiTcpisiordk  Ikace,  L^unov  Equation,  Dual  Norms. 

AMS(M08)  mOftet  etatifieation:  65F35. 

Abbreviatti  iMe  Lyapunov  Ikace  Bounds 
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Duu?e  thfraiiiT'r  ^ approximations  associated  with  com- 
bolt  tn!  function  are  shown  to  be  equivalent  to  iterations  involving  the  hyper- 

ohe  tengent  and  its  inverse.  Using  this  equivalent  formulation  many  results  about  these 
Fade  Iterations,  such  as  global  convergence,  the  semi-group  property  Sr  composi^^^^^ 

^  Obtained  easily.  In  the  second  part  of  the 
paper  it  is  shown  that  the  behavior  of  points  under  the  Fade  iterations  can  be  expressed  lia 

itera^on  Uesetwr;  “  t  of  »  completely  regular  iteration  and  a  chaotic 

itera  ion.  These  two  iterations  are  decoupled,  with  the  chaotic  iteration  taking  the  form  of  a 


1.  Introduction 

A  new  family  of  rational  iterations  was  introduced  in  [10]  for  computing  the 
matnx  sign  function.  These  iterations  are  based  on  Fade  approximatfons  for  a 
certain  hypergeometnc  func^  Table  1  of  [10]  gives  the  first  few  such  formulas 

noi  lu  H-ff numerator  and  denominator  poly- 
nomals  differ  m  degree  by  1  are  globally  convergent.  In  the  sequel  we  shall  refer 

Sml'Sp formulas  defined  by  such  rational  functions  as  the  principal  Fade 

ZTtlV  described  in  this  paper  will  be  carried  out 

for  the  scalar  ca^.  However,  it  must  be  emphasized  that  the  scalar  theory  carries 
over  directly  to  the  matrix  case  (see  [10]  for  details). 

ered“  ?!i'  “teresting  aspects  of  the  principal  Fade  iterations  are  consid- 

itT”  an  expression  of  these  iterations  in  terms  of  the  hyperbolic  tangent 
d  Its  inverse,  while  the  second  deals  with  chaotic  and  regular  behavior  of  points 

Foundation  under  Grant  No.  ECS- 
thi.  nffi  ’  f  7^*^  Scientific  Research  under  Grant  no.  F49620-94-1-0104DEF  and 
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Abstract 

A  survey  of  the  matrix  sign  function  is  presented,  including  the  historical  background, 
definitions  and  properties,  approximation  theory  and  computational  methods,  condition 
theory  and  estimation  procedures,  as  well  as  application  to  the  areas  of  control  theory, 
eigendecompositions  and  roots  of  matrices.  Many  new  results  are  presented  along  with  a 
description  of  the  Green’s  matrix  theory  developed  by  Godunov  and  others. 

1  Introduction 

The  matrix  sign  function  can  be  used  to  find  bases  for  the  positive  and  negative  invariant 
subspaces  of  a  matrix.  Since  many  important  problems  in  the  areas  of  control  theory  and 
eigenvalue-dgenvector  decompositions  have  solutions  that  can  be  given  in  terms  of  these  sub¬ 
spaces,  the  matrix  sign  function  has  enjoyed  renewed  interest  in  the  last  twenty  years.  This 
paper  summarizes  the  developments  of  this  period  and  then  looks  at  the  theory  and  applications 
of  the  sign  function. 

hi  1877  Zolotarjov  [98]  characterized  the  best  rational  approximations  of  the  scalar  sign 
function  in  terms  of  elliptic  sine  functions.  However,  this  result  does  not  seem  to  be  well  known 
and  the  introduction  of  the  matrix  sign  function,  especially  into  the  area  of  control  theory, 

*Tliii  icfeuch  was  supported  ia  part  by  the  Air  Force  Office  at  Sdeatific  Research  under  Grant  No.  F49620- 
94-1-0104DEF,  the  National  Science  Foundation  under  Grant  No.  ECS-9120S43  and  the  Office  of  Naval  Research 
under  Grant  No.  N00014-92-J-1706. 
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Abstract 

An  algorithm  is  presented  for  recovering  the  scaling  coefficients  from  a  given  scaling 
function.  This  procedure  is  then  extended  to  the  problem  of  finding  scaling  functions 
that  approximate  a  given  function. 


1  Introduction 

Given  Co,  Cl, ...,  Cn  such  that 

c**2,  (1) 

kmO 

define  the  associated  scaling  function  4  l>y  ihst  setting  4{^)  ^  0  for  x  outside  the  interval 
(0,n).  For  z  €  (0,n)  set 

^(*)  =  E  -  k) ,  (2) 

SaO 

subject  to  the  normalization  condition 

/  ^(z)dr  »  1  .  (3) 

/-OO 

Necessary  and  sufficient  conditions  for  ^  to  be  in  L^[0,n]  are  given  by  Daubechies  and 
Lagarias  [1]  and  a  irdl  known  dyadic  procedure  for  constructing  4{*)  et  rational  points 

2  €  (0, n)  from  the  values  ce,ei, . . . , c^  is  pven  by  Strang  in  [2].  See  aJbo  [3]. 

In  this  note  we  consider  the  inverse  problem:  if  we  are  given  a  scaling  function  4>  how 
do  we  recover  the  values  of  the  scaling  coefficients  cp,  ci, . . . ,  c^?  The  answer  to  this  ques¬ 
tion  opens  up  the  possibility  of  finding  scaling  functions  which  approximate  more  general 
functions. 

The  key  observation  is  that  co  can  be  recovered  byusing  (2)  for  z  in  the  interval  (0, 1/2). 
Since  ^(x)  s  0  for  x  <  0,  if  x  €  (0, 1/2)  then 

^(x)=cu^(2x).  (4) 

'This  research  was  supported  by  the  Office  of  Naval  Research  under  ONR  Grant  Number  NOOO 14-92- J- 
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